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The goal of the Advanced Formulation Development portion of the Inelastic Analysis Methods 
contract is the development of an alternative stress analysis tool, distinct from the finite 
element method, applicable to the engineering analysis of gas turbine engine structures. The 
boundary element method was selected for this development effort on the basis of its already 
demonstrated applicability to a variety of geometries and problem types characteristic of 
gas turbine engine components. This paper describes briefly major features of the BEST3D 
computer program and outlines some of the significant developments carried out as part the 
Inelastic Methods Contract. 


BEST3D OVERVIEW 

BEST3D (Boundary Element Stress Technology - Three Dimensional) is a general purpose 
three-dimensional structural analysis program utilizing the boundary element method. The 
method has been implemented for very general three-dimensional geometries, for elastic, in- 
elastic and dynamic stress analysis. Although the feasibility of many of the capabilities pro- 
vided had been demonstrated in a number of individual research efforts, the present code is 
the first in which they have been made available for large scale problems in a single code. 
In addition, important basic advances have been made in a number of areas, including the 
development and implementation of a variable stiffness plasticity algorithm, the incorporation 
of an embedded time algorithm for elastodymanics and the extensive application of particu- 
lar solutions within the boundary element method. Major features presently available in the 
BEST3D code include: 

• Very general geometry definition, including the use of doubly curved isoparametric sur- 
face elements and volume cells, with provision of full substructuring capability 

1 The work discussed in this paper was carried out as part of National Aeronautics and Space Administration 
contract NAS3-23697, “3-D Inelastic Analysis Methods for Hot Section Components”. The program manager 
at NASA-Lewis Research Center is Mr. C. C. Chamis. 
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• General capability for the definition of complex, time-dependent boundary conditions 

• Capability for nonlinear analysis using a variety of algorithms, solution proceedures and 
constitutive models 

• A very complete elastodynamic capability including provision for free vibration, forced 
response and transient analysis 

BEST3D has been successfully implemented on a variety of computers, including the IBM 
3090, various CRAY models and the Hewlett-Packard HP9000. 

The analytical basis and numerical implementation of the boundary element method for the 
major problem types considered are very briefly reviewed in the next two sections. Particu- 
lar attention is devoted to the variable stiffness plasticity algorithm and the time-embedded 
elastodynamic algorithm. Full details of both the analysis and implementation may be found 
in references 1, 2 and 3. 


QUASI-STATIC ANALYSIS 


By making use of the reciprocal work theorem, the governing differential equations for a three- 
dimensional (homogeneous) structure under combined thermal, mechanical and body force 
loadings can be converted to an integral equation written on the surface of the structure. 
This integral equation is: 

cijUi = [ (G tj u - F,j u , ) dS + [ (Gijfi + (3WjT) dV (1) 

J s Jv 

where T = temperature, W 3 — T tkj 6 x (3 = coefficient of thermal expansion, T tk j — the stress, 
a t k, due to a point force system, ej, and G XJ ,T tkj and F tJ are defined reference 1. The equation 

a t j = J(D t j k tk — S tjk u k )dS + Jj,T{j k t k + M XJ T)dV (2) 

allows calculation of stresses at any interior point where they are required. A similar equation 
for interior displacements can be obtained by setting c,y = 6 XJ in (1). 

In a purely elastic problem BEM stress analysis can be carried out entirely on the boundary 
of the structure. Once a physically reasonable set of boundary conditions has been prescribed, 
(1) can, in principle, be solved for all of the remaining boundary displacements and tractions. 

It is generally impossible to solve (1) exactly for real structures and loading conditions. Suit- 
able approximations of the boundary geometry, displacements and tractions must be used in 
order to reduce (1) to a system of algebraic equations. The present version of BEST3D mod- 
els boundary geometry and boundary values of field quantities using linear and/or quadratic 
isoparametric shape functions. The surface integrals in (1) are then evaluated numerically 
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using product Gaussian quadrature rules. The numerical implementation of the BEM is dis- 
cussed in detail in textbooks (ref. 4), as well as in references 1 and 2. 

In the case of inelastic analysis, the volume integrals in (1) cannot be calculated a priori , since 
they require knowledge of inelastic strain, which is itself a part of the solution. In this case 
equations (l), (2) and the inelastic material model can be regarded as a coupled system of 
nonlinear equations. In the numerical implementation of the BEM (2) is used to calculate the 
stresses at interior points, and the nonlinear material model is then used to evaluate inelastic 
strain. Since the volume integrals of inelastic strain vanish except in regions of nonlinear 
material response, approximations of geometry and field quantities are required only where 
nonlinearity is expected. In the original version of BEST3D, strain variation in the interior was 
represented using isoparametric volume cells, with the solution carried out using a relatively 
standard iteration proceedure. More recently, a new approach has been developed which 
exploits certain features of the constitutive relationships involved. The unknown nonlinear 
terms in the interior are now defined as scalar variables. A new direct numerical solution 
scheme comparable to the variable stiffness method used in finite element analysis has been 
developed and implemented, avoiding the requirement for an iterative solution. 


For a standard elasto-plastic flow problem the evolution of plastic flow is governed by: 


F{cr,j,h ) 

cf. 


A 


OF 

da. 


tj 


( 3 ) 

( 4 ) 


These equations together with the consistency relations (i.e., the stress point must remain 
on a newly developing yield surface characterized by a change in the hardening parameter h ) 
leads to an expression for the unknown plastic flow factor A as: 

a = r. (5) 

where 

J<r _ j_ dF 

'] H daij 

dF dh , dF 
dh de mn 

It should be noted that L*. depends upon the current state variable, not on the incremental 
quantities. 


mn 


However, the relationship given by (5) does not exist for ideal plasticity, as H vanishes for 
zero hardening. This can be avoided by reformulating the above expression in terms of strain 
increments: 

A = lUn (6) 

where 

J_ d£ 

dFj dF __ ( dF_ , 8FdF_. d£ 

da kl Vk,mn damn [ dc p kl + dh d( p k[ > da kl 



243 



where D tJ ki is the elastic constitutive tensor. It is evident that H' does not vanish for zero 
hardening (ideal plasticity). 

The basic boundary element formulation for an inelastic body undergoing infinitesimal strain 
is given by: 

= J r [GnM)U(x) - dT + J a B ijk (z,0 <r° dU (7) 

The stress rates at an interior point £ are obtained from equation (7) via the strain-displacement 
relations and the constitutive relationships = D*. u e k i — <r°.) as 

= j T [G’ ik (x,l,)u(x)-r it (z,f,)<H^)]d r (8) 

+ J n 0 ■ K '. + r <7^(0 

where the kernel functions have been defined in the reference 2. 

In equation (8) the volume integral must be evaluated in the sense of (Q — D ) with the limit 
taken as D — ► 0, where D is a spherical exclusion of small arbitrary radians with £ as its 
center. The term J a is the jump term derived from the analytical treatment of the integral 
over D . It is of considerable interest to note that the value of J * is independent of the size of 
the exclusion D , provided the initial stress distribution is locally homogeneous, i.e. uniform 
over its volume. 

The evaluation of strains and stresses at boundary points can be accomplished by considering 
the equilibrium of the boundary segment and utilizing constitutive and kinematic equations. 
The stresses and global derivatives of the displacements which lead to strains at a point £ can 
be obtained from the following set of coupled equations: 

- (A 6iju kt k{£) + /x(«ij(0 + Wj,, ■(£))) = (0 

= U(0 (9) 

duiij 1 _ duiiQ 

dm d£k d m 

where rji are a set of local axes at the field point £. 


All the above nonlinear formulations include initial stresses in the governing equations which 
are not known a priori and, therefore, are solved by using iterative procedures. A non-iterative 
direct solution procedure is made feasible in this work by reducing the number of unknowns 
in the governing equations by utilizing certain features of the incremental theory of plasticity 
expressed by equations (3) to (6). The initial stresses o^. appearing in equations (7) to (9) 
can be expressed in the context of an elastoplastic deformation as: 


&° = Kij\ 

ij J 


( 10 ) 


where Kij = D 


e dF_ 
bw da ki ‘ 
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Substituting (5) and (10) in equations (7) and (8) we can obtain: 


CijU,{£) = J[G t j(x)U(x) - F tJ (x,^)ui(x)]dT (11) 

and 

+1' ({) jf ^(*,0 K{ p (x) X(x) d(l (12) 

+ l' k J: pik K<M)k() 

Equations (11) and (12) can be solved simultaneously to evaluate the unknown values of 
displacements, traction rates and the scalar variable A. 

The equations for the boundary nodes (9) are similarly transformed to express them in terms 
of the scalar variable A using equations (5 ) and (10). 


TRANSIENT STRESS ANALYSIS 


The direct boundary integral formulation for a general, transient, elast.odynamic problem can 
be constructed by combining the fundamental point force solution of the governing equations 
(Stokes’ solution) with Grafh’s dynamic reciprocal theorem. Details of this construction can 
be found in Banerjee and Butterfield (ref. 4). For zero initial conditions and zero body forces, 
the boundary integral formulation for transient elastodynamics reduces to: 

*;(£M£.T) = J s [G tJ (x,£,T) * t t (x,T) - F t} (x,£,T) * mfaT)] dS(x) (13) 

where 

G,j*ti = f Gij(x,T]£, t) ti(x,r) dr (14) 

Jo 

fT 

F tJ * Ui = / Fij(x,T](, t) Ut(x,r) dr 

Jo ~ 

(15) 


are Reimann convolution integrals and £ and x are the space positions of the receiver (field 
point) and the source (source point). The fundamental solutions G{j and F XJ are the displace- 
ments and tractions at a point x and at a time T due to a unit force vector acting at a point 
£ at a time r. Equation (13) represents an exact formulation involving integration over the 
surface as well as the time history. It should also be noted that this is an implicit time-domain 
formulation because the response at time T is calculated by taking into account the history of 
surface tractions and displacements up to and including the time T. Furthermore, equation 
(13) is valid for both regular and unbounded domains. 
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Once the boundary solution is obtained, the stresses at the boundary nodes can be calculated 
without any integration by using the scheme described for the static case. For calculating 
displacements at interior points equation (13) can be used with c l} — 6, } and the interior 
stresses can be obtained from 


,(£,T) = (16) 


The functions G*. and F*. in the above equation are derived from the Stokes’ solution by 
differentiation. 


In the initial version of BEST3D, constant time stepping was used to obtain the transient 
dynamic response. It was found to be more effective to use a linear time variation of u and t 
on the boundaries. In this case: 


u t (x,T) = £[Mi u" \x) + M 2 u n t {x)} (IT) 

n= 1 

U(x,t) = ^2[Mi <" _1 (i,r) + M 2 <"($.)] (18) 

n — 1 

where M \ and M 2 are the time functions, and are of the form: 

Mi = ( 19 ) 

M 2 = ^=i^n(r) (20) 

For illustration purposes, consider the boundary integral equation for the first time step, i.e. 

CijUifoTi) - [ f [GijU - FijUi\ dS dr = 0 (21) 

Jt 0 Js 

The time integration in equation (21) by utilizing (18) is done analytically. After the usual 
numerical integration and assembly process, the resulting system equation is of the form: 

K][X‘] - - [B,'][r 0 ] = 0 ( 22 ) 


where: 

• A and B are matrices related to the unknown and known field quantities, respectively: 

• X and Y are the vectors of unknown and known field quantities, respectively: 

• for and Y the superscript denotes the time: 

• for A and B the superscript denotes the time step at which they are calculated, and the 
subscript denotes the local time nodes (1 or 2) during that time-stepping interval. 
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( 23 ) 


Since all the unknowns at time T = 0 are assumed to be zero, equation (22) reduces to: 

[A 1 2 )[X 1 ] = [B' 2 )[Y 1 ] + [B 1 l ][Y 0 } 

For second time step, the assembled system equation has the form 
KP' 2 ] - [Bl][Y 2 ] + [A\][X l ] - [B'JY 1 ] = -[^][P] + [BUY 1 ] - [A\][X°] + [B^Y 0 ] (24) 


As in the constant time variation scheme, only the matrices on the right hand side of equation 
(24) need be evaluated. However, one needs to integrate and assemble four matrices at each 
time step as compared to two in the case of constant time variation. This can be done with 
only a small increase in computing time by integrating all the kernels together and then 
assembling all the matrices together. Equation 24 can be rearranged such that: 

{A\]{X 7 ] = [B’][F J ] - [A\ + A'JX'} - [B\ + Bl)[Y' ] + [b|J[f‘] (25) 


In the above equation, all the quantities on the right hand side are known. Therefore, the 
unknown vector X at time T 2 can be obtained by solving the above equation. 

Thus, for the present case, the boundary integral equation (25) can be written in discretized 
form as: 


N 


[A\}[X N \ - [Bj][y"] = - ElK + - [b; + b"- , ][f"-'' + ') + [b"][f 0 ]] ( 26 ) 


n = 2 


or 

KP"] = [^r*] + [**] (27) 

It is of interest to note that, if time interpolation functions M 1 and M 2 are replaced by 
Mi pM 2 ~ 0.5</> n (r), the time stepping scheme for linear variation can be used for the case 
of constant variation with averaging between the local time nodes. 


SYMBOLS 


XL-t , ti 
Gij 

Fx 

fi 

T 

P 

k > Tj j k , Dij k ) \j k 

S 

V 


Kronecker delta symbol 
boundary displacements and tractions 
displacement point load solution 
traction kernel derived from G tJ 
mechanical body forces 
temperature 

coefficient of thermal expansion 
higher order kernels derived from G t j 
surface of three-dimensional structure 
interior of a three-dimensional structure 
stress tensor 
strain tensor 
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A 

F(a tJ ,h) 

h 


*, same as cr lJ 
superscript p , 
superscript e , 
D € ... 

ij kl 

x,y,z 

superscript °, 


* jump terms in boundary integral equation 
plastic flow factor 
yield function 
hardening parameter 
time derivative 
as in plastic component 

as in elastic component 

elastic constitutive tensor 
points in three-dimensional space 
as in <7°. initial stress (or strain) 




higher order kernel derived from G tJ 


A 

i,r 


*, as in G X] * t t 

M\ , M2 


[At], etc. 


Lame constants 

denote time in dynamic analysis 
time convolution 

shape functions for time variation 
coefficient matrices in time domain solution 
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